NAVAL 

POSTGRADUATE 

SCHOOL 

MONTEREY,  CALIEORNIA 


THESIS 


PERFORMANCE  OF  THE  SELF  REFERENCING 
INTERFEROMETER  IN  THE  PRESENCE  OF 
SIMULATED  DEEP  TURBULENCE  AND  NOISE 
EFFECTS 

by 

Lee  T.  Johnson 
December  2013 

Thesis  Co-Advisors:  Brij  N.  Agrawal 

Jae-Jun  Kim 


Approved  for  public  release;  distribution  is  unlimited 


THIS  PAGE  INTENTIONALLY  LEET  BLANK 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 
Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instruction, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send 
comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to 
Washington  headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA 
22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188)  Washington  DC  20503. 

2.  REPORT  DATE 

December  20 1 3 


6.  AUTHOR(S)  Lee  T.  Johnson 


11.  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official  policy 
or  position  of  the  Department  of  Defense  or  the  U.S.  government.  IRB  protocol  number _ N/A 


13.  ABSTRACT  (maximum  200  words) 

Current  laser  weapon  systems  are  limited  to  close  range  encounters  because  the  laser  beam  attenuates  quickly  within 
the  atmosphere.  A  phenomenon  known  as  deep  turbulence  is  characterized  by  strong  scintillation  and  branch  points 
in  the  wave-front  phase.  Many  wave-front  sensors  perform  poorly  in  the  presence  of  deep  turbulence,  and  are  unable 
to  accurately  reconstruct  the  wave-front.  This  paper  examines  a  wave-front  sensor,  the  self-referencing  interferometer 
(SRI)  that  is  theoretically  immune  to  the  effects  of  deep  turbulence. 

The  SRI  is  both  simulated  mathematically  and  constructed  in  the  lab  for  comparison  between  analytical  and 
experimental  results.  Performance  of  the  SRI  is  analyzed  in  the  presence  of  realistic  deep  turbulence  effects 
generated  by  a  spatial  light  modulator,  and  realistic  noise  effects  introduced  by  the  digital  imaging  system.  Simulated 
results  show  a  significant  loss  of  signal  level  as  turbulence  is  increased,  but  a  resilience  of  the  wave-front  sensor 
above  a  signal-to-noise  ratio  of  two.  Analogously,  in  the  experimental  results  the  signal  drops  off  rapidly  with 
increasing  levels  of  turbulence,  and  reaches  unacceptably  low  levels  above  a  Rytov  number  of  0.4.  A  qualitative 
analysis  of  the  wave-front  reconstruction  shows  remarkable  similarity  between  simulated  and  experimental  results, 
though  the  experimental  results  contain  far  more  error  induced  branch  points  than  in  the  simulation.  Methods  are 
being  explored  to  boost  the  signal  and  reduce  the  noise  at  the  camera  to  allow  the  system  to  handle  higher  levels  of 
turbulence. 


16.  PRICE  CODE 


NSN  7540-01-280-5500  Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  239-18 


20.  LIMITATION  OF 
ABSTRACT 


15.  NUMBER  OF 
PAGES 

67 


14.  SUBJECT  TERMS  Adaptive  optics,  self-referencing  interferometer,  spatial  light  modulator,  deep 
turbulence 


18.  SECURITY 
CLASSIFICATION  OF  THIS 
PAGE 

Unclassified 


19.  SECURITY 
CLASSIFICATION  OF 
ABSTRACT 

Unclassified 


17.  SECURITY 
CLASSIFICATION  OF 
REPORT 

Unclassified 


12b.  DISTRIBUTION  CODE 

A 


12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Postgraduate  School 
Monterey,  CA  93943-5000 

9.  SPONSORING  /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

N/A 


5.  FUNDING  NUMBERS 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


4.  TITLE  AND  SUBTITLE 

PERFORMANCE  OF  THE  SELF  REFERENCING  INTERFEROMETER  IN  THE 
PRESENCE  OF  SIMULATED  DEEP  TURBULENCE  AND  NOISE  EFFECTS 


3.  REPORT  TYPE  AND  DATES  COVERED 

Master’s  Thesis 


1.  AGENCY  USE  ONLY  (Leave  blank) 


1 


THIS  PAGE  INTENTIONALLY  LEET  BLANK 


11 


Approved  for  public  release;  distribution  is  unlimited 


PERFORMANCE  OF  THE  SELF  REFERENCING  INTERFEROMETER  IN  THE 
PRESENCE  OF  SIMULATED  DEEP  TURBULENCE  AND  NOISE  EFFECTS 


Lee  T  Johnson 

Lieutenant,  United  States  Navy 
B.S.,  North  Carolina  State  University,  2009 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  ASTRONAUTICAL  ENGINEERING 

from  the 


NAVAL  POSTGRADUATE  SCHOOL 
December  2013 


Author: 


Lee  T.  Johnson 


Approved  by:  Brij  N.  Agrawal 

Thesis  Co-Advisor 


Jae-Jun  Kim 
Thesis  Co-Advisor 


Knox  T.  Millsaps 

Chair,  Department  of  Meehanieal  and  Aerospaee  Engineering 


THIS  PAGE  INTENTIONALLY  LEET  BLANK 


IV 


ABSTRACT 


Current  laser  weapon  systems  are  limited  to  elose  range  encounters  because  the  laser 
beam  attenuates  quickly  within  the  atmosphere.  A  phenomenon  known  as  deep 
turbulence  is  characterized  by  strong  scintillation  and  branch  points  in  the  wave-front 
phase.  Many  wave-front  sensors  perform  poorly  in  the  presence  of  deep  turbulence,  and 
are  unable  to  accurately  reconstruct  the  wave-front.  This  paper  examines  a  wave-front 
sensor,  the  self-referencing  interferometer  (SRI)  that  is  theoretically  immune  to  the 
effects  of  deep  turbulence. 

The  SRI  is  both  simulated  mathematically  and  constructed  in  the  lab  for 
comparison  between  analytical  and  experimental  results.  Performance  of  the  SRI  is 
analyzed  in  the  presence  of  realistic  deep  turbulence  effects  generated  by  a  spatial  light 
modulator,  and  realistic  noise  effects  introduced  by  the  digital  imaging  system. 
Simulated  results  show  a  significant  loss  of  signal  level  as  turbulence  is  increased,  but  a 
resilience  of  the  wave-front  sensor  above  a  signal-to-noise  ratio  of  two.  Analogously,  in 
the  experimental  results  the  signal  drops  off  rapidly  with  increasing  levels  of  turbulence, 
and  reaches  unacceptably  low  levels  above  a  Rytov  number  of  0.4.  A  qualitative  analysis 
of  the  wave-front  reconstruction  shows  remarkable  similarity  between  simulated  and 
experimental  results,  though  the  experimental  results  contain  far  more  error  induced 
branch  points  than  in  the  simulation.  Methods  are  being  explored  to  boost  the  signal  and 
reduce  the  noise  at  the  camera  to  allow  the  system  to  handle  higher  levels  of  turbulence. 
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I.  INTRODUCTION 


A.  BACKGROUND  AND  MOTIVATION 

High  energy  laser  (HEL)  weapon  systems  have  been  a  focus  of  research  for  over 
25  years  in  the  Department  of  Defense.  Laser  weapon  systems  have  historically  been 
researched  by  the  Air  Force  as  a  solution  for  missile  defense,  and  have  been  tested  aboard 
aircraft.  Though  tested  successfully  to  defend  against  missiles  at  short  range  in  the  early 
1980s  aboard  the  Airborne  Laser  Laboratory  (ALL),  the  system  was  never  introduced 
operationally.  Many  advances  in  the  field  of  beam  control  and  adaptive  optics  have 
generated  interest  in  a  longer-range  laser  weapon  system  to  shoot  down  ballistic  missiles 
in  boost  phase.  The  Airborne  Laser  Test  Bed  (ALTB)  program  was  established  in  the 
mid-1990s  to  explore  this  possibility.  That  program  was  cancelled  in  2011  after  more 
than  a  decade  of  research  and  development,  but  no  operational  system  is  yet  being  fielded 
as  a  result  of  the  effort  [1],  [2].  In  recent  years  the  U.S.  Navy  has  been  actively  pursuing 
the  development  of  laser  weapon  system  for  use  on  ships.  There  have  been  some 
successes  in  the  Navy’s  laser  weapon  system  (LaWS)  research,  and  the  system  has  been 
employed  to  disable  some  small  unmanned  aircraft  and  small  boats.  That  system  is  now 
scheduled  to  be  installed  aboard  a  navy  ship  in  2014  for  operational  use  as  a  technology 
demonstrator  [3],  [4].  New  laser  technologies  are  advancing,  and  will  be  considered  for 
various  military  applications  in  the  next  decade.  These  new  technologies  include  solid- 
state  lasers,  and  free-electron  lasers.  New  applications  for  HEL  weapon  systems  may 
include  air  defense,  ground  and  maritime  offense,  space  control,  and  urban  operations  [5]. 

While  there  are  differences  in  each  of  these  laser  applications  and  technologies, 
there  are  common  elements.  Military  applications  of  lasers  would  benefit  greatly  from 
increased  range  and  accuracy,  whether  the  laser  system  is  in  the  kilowatt  or  megawatt 
range.  All  laser  beams  suffer  from  atmospheric  attenuation  when  propagating  through 
the  air.  This  attenuation  is  especially  troublesome  in  the  maritime  environment  because 
of  increased  particulates  in  the  air  and  relatively  horizontal  path.  While  the  particles 
within  the  air  and  air  molecules  themselves  can  absorb,  reflect,  scatter,  and  otherwise 
disperse  the  energy  contained  in  the  laser  beam,  the  motion  of  the  air  and  the  thermal 
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gradients  between  air  masses  will  also  eause  a  beam  to  lose  its  eoherenee  as  it  passes 
through  this  turbulent  medium.  This  is  eaused  by  rapidly  ehanging  indiees  of  refraetion, 
and  inhomogeneities  in  the  small  masses  of  air  that  the  beam  eneounters  along  its  path 
[6].  A  defining  eharaeteristie  of  laser  beams  is  their  spatial  eoherenee.  Spatial  eoherenee 
is  eharaeterized  by  a  eonstant  wave-front  amplitude  and  phase  profile,  whieh  allows  a 
laser  beam  to  maintain  its  small  diameter  over  long  distances.  If  a  laser  beam  encounters 
a  turbulent  medium  it  will  lose  its  coherence  [7].  A  loss  of  coherence  in  the  beam  will 
cause  beam  dispersion,  and  the  energy  contained  in  the  beam  will  be  spread  out.  The 
farther  the  beam  travels  through  the  turbulent  medium,  the  more  aberrated  the  beam  will 
become.  The  turbulent  effects  will  eventually  produce  scintillation  in  the  beam,  which  is 
an  amplitude  fluctuation  across  the  electromagnetic  field.  This  scintillation  causes  rapid 
degeneration  of  the  beam  and  very  little  energy  makes  it  to  the  target.  The  field  of 
adaptive  optics  is  dedicated  to  the  study  of  correcting  aberrated  wave-fronts. 

Many  developments  have  been  made  in  correcting  for  weak  turbulence. 
Astronomical  seeing  from  ground  based  telescopes  is  now  very  near  the  diffraction  limit 
due  to  the  advancement  of  adaptive  optics  [8].  While  these  advances  have  greatly 
impacted  the  scientific  community,  there  is  still  much  to  be  done  in  the  area  of  HEL 
weapon  systems.  Weapon  systems  that  would  often  benefit  from  increased  range  and 
accuracy  within  the  Earth’s  atmosphere  have  a  much  bigger  challenge  than  the  relatively 
short  path  lengths  and  weak  turbulence  encountered  by  ground-based  skyward-  looking 
telescopes.  This  research  focuses  on  adaptive  optics  applications  requiring  long  path 
lengths  through  a  turbulent  medium. 

B,  DEEP  TURBULENCE 

Eor  the  purposes  of  this  paper,  the  term  “deep  turbulence”  will  be  used  to  describe 
aberrations  with  significant  scintillation  and  branch  points  in  the  wave-front.  Adaptive 
optics  is  still  advancing  in  the  area  of  deep  turbulence.  Traditional  wave-front  sensors 
such  as  the  lateral  sheering  interferometer  and  Shack-Hartmann  wave-front  sensor  are 
highly  sensitive  to  the  presence  of  scintillation  and  branch  points,  but  there  are  more 
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novel  wave-front  sensors  that  are  naturally  immune  to  the  effeets  of  seintillation  and 
branch  points  [9].  One  such  sensor,  the  SRI,  is  the  subject  of  this  investigation. 

1.  Scintillation 

Scintillation  is  the  fluctuation  in  amplitude  of  an  optical  wave-front  caused  by 
propagation  through  long  path  lengths  or  strong  turbulence  (or  both)  within  a  turbulent 
medium.  When  scintillation  is  strong  enough,  there  can  be  areas  of  very  low  to  zero 
intensity  (squared  amplitude)  in  the  wave-front  that  present  challenges  to  wave-front 
reconstructors.  Scintillation  is  caused  by  the  constructive  and  destructive  interference  of 
individual  rays  of  light  from  a  point  source  as  they  pass  through  a  medium  with  varying 
indices  of  refraction  between  individual  pockets  along  the  propagation  path.  These  minor 
angular  deflections  of  individual  light  rays  worsen  over  long  path  lengths  to  produce 
significant  aberrations  in  the  wave-front. 

2,  Branch  Points  and  Branch  Cuts 

When  scintillation  becomes  severe,  as  is  the  case  with  deep  turbulence,  the 
amplitude  function  can  drop  low  enough  to  create  areas  of  zero  intensity  in  the  wave- 
front.  These  areas  represent  a  discontinuity  in  the  phase  function  of  the  reconstructed 
wave-front.  The  discontinuity  is  physically  interpreted  as  a  jump  of  27i  radians  in  the 
phase.  These  discontinuities  are  known  as  branch  points,  and  are  unavoidable  with 
increasing  levels  of  turbulence  [10].  Branch  points  occur  in  oppositely  signed  pairs,  and 
are  connected  by  an  arbitrarily  positioned  line  called  a  branch  cut.  Branch  points  are 
generally  linked  by  a  trough  of  low  intensity.  This  trough  is  the  most  natural  place  to 
position  the  branch  cut.  This  27i  jump  causes  a  major  problem  for  many  types  of  wave- 
front  reconstructors,  and  continuous  face-sheet  deformable  mirrors,  which  have 
unavoidable  correction  error  in  the  vicinity  of  branch  cuts.  Beam  correction  is  most 
important  in  areas  of  high  intensity,  and  thus  branch  cuts  are  placed  in  areas  of  lowest 
intensity  [10].  Traditional  wave-front  reconstructors,  such  as  the  least  mean  square  error 
wave-front  reconstructor,  calculate  the  phase  of  a  point  of  interest  by  summing  the  phase 
differences  along  all  possible  paths  from  a  phase  reference  point  to  the  point  of  interest, 
and  then  averaging  the  multiplicity  of  results  to  reduce  the  effects  of  noise.  This  method 
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works  well  unless  a  number  of  those  paths  eross  a  braneh  eut,  and  a  2n  diseontinuity 
skews  the  results  of  the  averaging  [11], 


3,  Rytov  Number 


The  Rytov  number  is  a  single  number  that  ean  be  used  to  describe  the  severity  of 
aberration  experienced  by  a  monochromatic  beam  as  it  propagates  through  a  turbulent 
medium.  This  number  is  useful  because  it  gives  an  idea  of  the  characteristics  of  the 
wave-front  one  should  expect.  Generally  a  Rytov  number  around  0.25  marks  the  onset  of 
strong  scintillation,  and  the  appearance  of  branch  points  in  the  wave-front  [12].  For  the 
purposes  of  this  paper,  aberrations  with  a  Rytov  number  above  0.25  will  be  considered 
deep  turbulence.  Rytov  number  is  a  function  of  propagation  path  length  (L),  wavelength 

of  light  (k),  and  the  index  of  refraction  structure  constant  which  is  a  measure  of  the 

severity  of  turbulence  in  the  medium.  The  Rytov  number  is  expressed  as  a  statistical 
variance  defined  as  [13] 


1.23C 


flTrX^  - 

—  £6 

y  A  J 


(1) 


Propagation  path  length  generally  varies  between  10s  and  100s  of  kilometers. 
Severe  turbulence  generally  has  a  value  greater  thanlO If  propagation 


path  length  is  short  enough  this  may  not  cause  scintillation  in  the  wave-front.  As  a 
corollary,  even  if  is  relatively  low  and  propagation  path  is  very  long,  scintillation  may 
still  appear.  This  is  why  Rytov  has  been  selected  as  a  measure  of  aberration  rather  than 
alone. 


C.  SELF-REFERENCING  INTERFEROMETER 

Classifying  the  nature  of  the  aberration  on  an  incoming  wave-front  involves  a 
camera,  which  acts  as  the  wave-front  sensor,  and  a  reconstructor,  which  interprets  the 
data  collected  by  the  wave-front  sensor  to  reconstruct  the  wave-front  phase.  This  paper 
analyzes  one  particular  kind  of  wave-front  sensor  called  the  self-referencing 

interferometer  (SRI).  The  SRI  is  ideal  for  use  with  deep  turbulence  because  it  measures 
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the  wave-front  field  direetly,  and  therefore  is  unaffected  by  branch  points  and 
scintillation.  The  SRI  accomplishes  direct  measurement  of  the  phase  in  a  unique  way 
that  does  not  require  the  use  of  a  traditional  reconstruction  method  such  as  a  least  mean 
square  error  reconstructor. 


The  self-referencing  interferometer  is  composed  of  a  series  of  optical  elements 
that  generate  four  phase-shifted  signal  beams  that  have  each  been  interfered  with  an 
equally  phase-shifted  reference  beam.  These  four  individual  beams  are  then  collected 
with  a  digital  imaging  system.  The  reconstruction  is  then  accomplished  with  a  simple 
formula  that  depends  only  on  the  intensity  profile  of  the  light  in  each  of  the  four  beams. 
Equation  2  is  used  to  calculate  the  phase  at  each  point  (/,  j)  in  the  wave-front  using  the 
four  individual  intensity  frames.  There  are  multiple  methods  for  creating  an  SRI  using 
more  or  less  than  four  frames.  This  particular  method  is  known  as  the  4-bucket  method. 
The  wave-front  phase  is  calculated  by  [14] 


(2) 


Where  h,  h,  h,  and  h  correspond  to  intensities  with  0,  nil,  k,  and  3k/2  phase  shift, 
respectively.  The  phase  output  will  be  from  -n  to  n  and  will  need  to  be  unwrapped  before 
using  a  continuous  face-sheet  DM  [15]. 


D.  THESIS  OBJECTIVE 

The  objective  of  this  research  is  to  create  and  characterize  a  SRI  at  the  Naval 
Postgraduate  School  for  use  in  on-going  adaptive  optics  and  beam  control  research.  The 
SRI  will  be  simulated  in  the  MATLAB  software  program,  and  constructed  in  the 
laboratory.  The  SRI  will  be  used  to  reconstruct  wave-fronts  containing  strong 
scintillation  and  branch  points.  The  performance  of  the  SRI  will  be  analyzed,  and  its 
limitations  probed.  Specifically,  the  ability  of  the  SRI  to  qualitatively  reproduce  wave- 
fronts  and  identify  branch  point  locations  will  be  assessed.  The  SRI  will  be  tested  at  a 
variety  of  turbulence  and  noise  levels  to  determine  the  limits  of  its  operation.  Finally, 
areas  will  be  identified  for  the  improvement  of  performance  of  the  SRI. 
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II.  SIMULATION 


A.  SELF  REFERENCING  INTERFEROMETER  SIMULATION 

The  self-referencing  interferometer  was  simulated  in  MATLAB  using  some  basic 
mathematical  equations  and  raw  data  provided  by  WaveProp.  WaveProp  is  MATLAB 
third  party  toolbox  written  by  the  Optical  Sciences  Company.  It  is  used  to  accurately 
simulate  the  effects  of  atmospheric  turbulence  on  passing  optical  wave-fronts.  The 
purpose  of  the  simulation  is  to  establish  a  model  against  which  the  experimental  set-up 
can  be  measured.  The  simulation  uses  many  of  the  same  processes  developed  for  the 
analysis  of  experimental  data,  and  addresses  the  reality  gap  by  implementing  realistic 
camera  noise.  The  simulation  provides  a  standard  of  performance  that  should  be 
attainable  in  experimentation  after  refinement  of  optical  alignment,  removal  of  excess 
system  noise,  and  optimization  of  data  processing  techniques.  Error  can  enter  into  the 
reconstruction  during  data  processing  because  of  a  phenomenon  known  as 
misregistration.  This  occurs  when  there  is  misalignment  between  the  four  frames  when 
the  wave-front  phase  is  calculated.  For  the  reconstruction  to  be  accurate,  each  pixel  of 
each  frame  must  align  with  its  counterpart  in  the  adjacent  frame,  so  that  the  phase  at  that 
location  can  be  calculated  from  the  correct  pixels. 

WaveProp  is  provided  a  value,  propagation  path  length,  and  wavelength.  It 

produces  amplitude  and  phase  data  corresponding  to  a  specific  Rytov  number.  Four 
intensity  profiles  are  then  generated  via  a  simple  formula  that  simulates  the  interference 
and  phase  shifting  of  the  given  amplitude  and  phase  data  [16].  These  frames  can  be 
calculated  using  the  equations, 

/j  =  A  +  A  cos  (p 
1^=  A-  Asincp 
=  A- A  cos  (p 
I^  =  A  +  Asmcp 

Once  the  four  frames  have  been  generated  the  phase  reconstruction  can  be  completed 
using  Equation  2. 
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B, 


CAMERA  NOISE 


The  simulated  eamera  noise  is  composed  of  several  different  types  of  noise.  Dark 
current  is  simulated  by  adding  a  frame  of  uniform  intensity  that  corresponds  to  an 
integration  time  and  a  fixed  rate  of  accumulation  of  electrons  based  on  sensor 
temperature.  Dark  current  can  be  very  small  (on  the  order  of  1  electron  per  second)  for  a 
good  imaging  system.  Imaging  systems  with  built  in  cooling  can  achieve  50°  C  below 
ambient  combined  with  integration  times  in  the  microsecond  range.  Pixel  non-uniformity 
is  simulated  by  generating  a  frame  of  Gaussian  distributed  numbers  with  a  standard 
deviation  of  some  small  percent,  and  a  mean  value  of  one.  This  frame  is  then  used  as  a 
scaling  factor  to  multiply  by  the  image  frame  to  simulate  the  variation  in  sensitivity  from 
pixel  to  pixel.  Both  dark  current  noise  and  pixel  non-uniformity  have  been  fixed  for  this 
simulation,  and  will  not  vary  with  changing  SNR.  For  this  simulation  a  dark  current  value 
of  10  electrons  per  pixel  per  second  was  used  at  an  integration  time  of  1x10'  seconds, 
and  the  value  for  pixel  non-uniformity  was  five  percent  standard  deviation. 

The  system  read  noise,  or  on-chip  noise,  represents  the  noise  introduced  by  the 
charge  coupled  device  (CCD).  This  value  can  be  very  low  for  a  low  noise  camera,  on  the 
order  of  15  photo  electrons  RMS  [17].  Read  noise  is  given  by  the  equation, 

=(««)'  (4) 

The  off-chip  noise  is  introduced  by  the  analog  to  digital  converter  (ADC),  and 
corresponds  to  the  error  introduced  in  the  conversion  of  photons  to  digital  signal  [18]. 
The  variance  in  this  noise  is  given  by  the  equation, 

(5) 

where  LSB  is  the  number  of  electrons  per  least  significant  bit  of  the  16-bit  ADC. 

Finally,  the  shot  noise  is  a  function  of  the  random  arrival  time  of  photons  within 
the  imaged  wave-front  and  is  described  by  the  equation, 
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where  is  the  quantum  efficieney  of  the  camera,  Nsput  is  the  number  of  interfero grams 

(four  in  this  case),  and  h  and  h  are  the  average  intensities  of  the  signal  and  reference 
beams  respectively  [18].  In  practice  the  signal  and  reference  beams  can  be  imaged 
separately  by  blocking  one  or  the  other  beam  to  get  an  accurate  value  for  the  intensities. 


This  value  will  obviously  change  depending  on  the  level  of  turbulence  [18].  The  total 
system  noise  is  then  given  by  equation. 


(7) 


This  standard  deviation  of  the  total  system  noise  can  be  used  to  generate  a  noise  frame 
Gaussian  noise  frame  is  then  added  to  the  image  frame.  These  three  types  of  noise  are 
used  to  realistically  simulate  the  noise  introduced  by  the  imaging  system,  and  can  be 
individually  tweaked  to  yield  a  signal  to  noise  ratio  (SNR)  comparable  to  the 
experimental  setup,  or  to  test  the  performance  under  different  noise  conditions. 

The  SNR  is  defined  as  the  mean  signal  strength  (//  )  divided  by  the  standard 
deviation  of  the  noise  ( ),  and  is  given  by  the  equation. 


SNR  =  -^ 


(8) 


While  this  is  obviously  only  an  approximation  of  a  much  more  complicated  and  dynamic 
effect,  it  can  be  used  to  measure  the  response  of  the  SRI  in  low  SNR  environments. 
Shown  in  Figure  1  is  a  sample  noise  frame  representative  of  noise  that  will  be  added  to 
test  the  SRFs  response  to  low  SNR. 
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Figure  1 .  Sample  noise  frame  used  to  test  SRI  resilience  to  Gaussian  noise  effects. 
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c. 


PERFORMANCE  METRIC 


The  quality  of  the  reeonstruction  is  measured  by  the  estimation  Strehl  ratio  [18], 


-l2 


(9) 


If  the  wave-front  is  reeonstructed  exaetly  as  it  was  entered,  the  Strehl  ratio  will  be 
one.  With  simulated  noise  effects  added  to  the  data,  the  Strehl  ratio  will  be  less  than  one. 
By  testing  a  variety  of  cases,  the  performance  of  the  SRI  can  be  estimated  in  the  presence 
of  varying  levels  of  turbulence  and  noise  effects.  The  original  wave-front  ( )  that  was 


generated  by  WaveProp  is  provided  to  the  SRI  simulation,  and  is  used  as  a  baseline 
against  which  the  reconstructed  wave-front  ( )  is  measured.  The  WaveProp  amplitude 


and  phase  data  are  scaled  to  match  the  frame  size  of  the  experimental  data  at  200  x  200 
pixels,  and  every  pixel  in  that  frame  has  its  phase  calculated  within  the  SRI  simulation  to 
yield  the  reconstructed  wave-front  data. 


The  simulation  was  run  for  several  values  of  Rytov  numbers  to  observe  the 
general  trends  of  increasing  both  Rytov  number  and  SNR  as  shown  in  Table  1. 


Rytov  Number 

SNR 

Strehl  Ratio 

0.301 

4 

0.950 

0.301 

2 

0.827 

0.301 

1 

0.539 

2.092 

4 

0.943 

2.092 

2 

0.822 

2.092 

1 

0.570 

3.313 

4 

0.942 

3.313 

2 

0.822 

3.313 

1 

0.568 

Table  1 .  Strehl  ratio  shown  for  a  variety  of  Rytov  numbers  and  SNRs. 
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The  trend  here  indicates  that  SRI  performance  is  theoretically  independent  of 
turbulence  level,  but  is  more  explicitly  linked  to  SNR.  It  is  true  that  in  most  cases 
increased  levels  of  turbulence  will  in  fact  attenuate  the  signal  and  lead  to  lower  SNR,  but 
if  signal  level  can  be  maintained  despite  high  levels  of  turbulence,  then  the  reconstruction 
retains  fidelity.  Also,  it  can  be  seen  that  for  Gaussian  noise  only,  the  SNR  cannot  drop 
much  below  two  for  the  reconstruction  to  remain  valid. 
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III.  EXPERIMENTATION 


A,  EXPERIMENTAL  SET-UP  AND  DATA  COLLECTION 

The  SRI  was  built  within  the  Adaptive  Opties  Center  of  Exeellence  (AOCoE)  at 
the  Naval  Postgraduate  School  (NPS).  The  major  components  of  the  setup  are  the  laser 
source,  the  spatial  light  modulator  (SLM),  the  SRI,  and  the  digital  imaging  system.  The 
experimental  setup  begins  at  the  source  of  the  laser  light.  In  this  case  the  source  is  a  2.5 
mW  laser  at  632  nm  (red).  The  laser  propagates  from  the  source  to  the  SEM  where  the 
wave-front  is  inflicted  with  a  high  degree  of  aberration  to  simulate  the  effects  of  deep 
turbulence.  The  reflected  light  from  the  SLM  is  polarized  so  that  it  follows  a  different 
return  path  than  the  incident  light.  The  light  is  then  propagated  to  the  SRI  where  it  is 
split  into  the  four  phase-shifted  interfered  beams  used  for  wave-front  reconstruction. 
These  components  are  part  of  a  more  complex  set-up  at  the  Adaptive  Optics  Center  of 
Excellence  (AOCoE)  at  the  Naval  Postgraduate  School  (NPS),  shown  in  Eigure  2. 


Eigure  2.  Adaptive  optics  test-bed  at  the  Naval  Postgraduate  School.  Optical  path  used 

in  SRI  experimentation  is  shown  here  in  red. 
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Figure  3  is  a  photograph  of  the  test-bed  as  it  is  today.  This  image  was  captured 
from  the  right  side  of  the  test-bed,  with  the  laser  source  just  off  camera  to  the  right. 


Figure  3.  Adaptive  optics  test-bed  at  the  Naval  Postgraduate  School. 


Increasing  the  level  of  turbulence  at  the  SLM  causes  the  signal  to  attenuate 
proportionally,  and  thus  there  is  a  loss  of  intensity  at  the  wave-front  sensor  as  Rytov 
number  increases.  This  loss  of  signal  causes  a  decrease  in  signal-to-noise  ratio,  which 
negatively  affects  the  overall  fidelity  of  the  reconstruction.  The  goal  is  to  either  reduce 
the  noise  enough,  or  increase  the  signal  enough  to  avoid  significant  impact  to  the 
reconstruction  process.  The  signal  strength  can  be  increased  by  replacing  low  quality 
optics  with  higher  quality,  more  reflective  optics  or  by  increasing  the  integration  time  at 
the  sensor.  The  noise  can  be  reduced  either  by  reducing  the  temperature  of  the  sensor, 
filtering  or  blocking  stray  light,  or  by  using  more  advanced  techniques  such  as  fiat 
fielding  or  frame  averaging  to  reduce  the  overall  level  of  Gaussian  random  noise  [19]. 

B.  SPATIAL  LIGHT  MODULATOR 

The  SLM  is  a  512x512  liquid  crystal  modulator  backed  by  a  dielectric  mirror. 
The  specific  model  used  is  the  Boulder  Non-Linear  Systems  XY  Ferroelectric  Series  with 
amplitude  only  modulation.  For  the  purposes  of  this  research,  the  SLM  was  used  to 
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generate  static  turbulence  for  characterization  of  the  SRI,  but  it  also  has  potential  for  use 
in  closed-loop  adaptive  optics.  The  SLM  is  capable  of  up  to  90  degrees  of  polarization 
rotation  at  1  kHz  frame  rate.  An  image  of  the  SLM  is  shown  in  Figure  4. 


Figure  4.  Spatial  light  modulator  manufactured  by  Boulder  Non-Linear  Systems  used  in 

the  experimental  set-up. 


WaveProp  was  used  to  generate  amplitude  and  phase  data  for  simulated  deep 
turbulence  effects,  including  scintillation  and  branch  points.  This  amplitude  and  phase 
data  is  then  implemented  by  the  spatial  light  modulator  (SLM)  to  aberrate  the  incident 
laser  light  [20].  The  light  reflected  from  the  SLM  is  then  representative  of  light  that  has 
passed  through  a  long  propagation  path  within  a  turbulent  medium.  For  the  experimental 
data  collection,  several  arbitrarily  selected  Rytov  numbers  were  used  ranging  from  0.044 
to  0.410.  The  intent  was  to  analyze  the  SRI  using  a  wide  range  of  values  to  determine  if 
there  were  any  trends  in  the  behavior  of  the  SRI  at  different  turbulence  or  signal  levels. 
The  performance  of  the  SLM  in  reproducing  the  turbulence  correctly  was  not  explicitly 
measured,  but  rather  was  taken  as  a  potential  source  of  error  in  the  final  reconstruction  at 
the  SRI. 

C.  SELF-REFERENCING  INTERFEROMETER 

Shown  in  Figure  5  is  a  detailed  view  of  the  SRI  setup  [21]. 
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Figure  5.  Optical  layout  of  the  self-referencing  interferometer  in  the  AOCoE  at  NFS. 


A  full  description  of  this  instrument  was  presented  at  the  2013  DBFS  Beam 
Control  Conference  in  Monterey,  California  [21].  The  SRI  is  a  polarization  phase 
shifting  interferometer  based  on  a  Mach-Zehnder  interferometer  [22],  [23].  The 
incoming  light  is  divided  by  a  polarized  beam  splitter  into  two  beams,  a  reference  beam 
and  a  test  beam.  The  reference  beam  is  generated  by  placing  a  pinhole  on  the  beam  path 
that  reflects  from  the  FBSl  to  spatially  filter  the  incoming  light.  Lens  LI  is  used  to 
collimate  the  reference  beam.  Lenses  L2  and  L3  are  used  to  relay  the  pupil  plane  of  the 
Test  Beam.  The  beams  pass  through  a  50/50  beam  splitter  and  the  two  output  beams,  A 
and  B,  pass  through  quarter  wave-plates.  Lenses  L4,  L6  and  L5,  L7  are  used  to  relay  the 
pupil  plane  of  beams  A  and  B,  respectively.  Both  beam  A  and  beam  B  are  composed  of  a 
reference  beam  and  a  test  beam.  Beam  A  passes  through  two  quarter  wave  plates  Q1  and 
Q2,  while  beam  B  passes  through  a  single  quarter  wave-plate,  Q3.  The  fast  axis  of  quarter 
wave-plate  Q1  is  set  at  0°  and  adds  a  %I2  phase  shift.  The  fast  axis  of  quarter  wave-plates 
Q2  and  Q3  are  set  at  45°.  Wave-plate  Q1  converts  reference  and  test  beam  of  beam  A  in 
to  opposite  handed  circularly  polarized  light.  Wave-plate  Q3  converts  beam  B  in  to 
opposite  handed  circularly  polarized  light.  The  additional  FBSl  splits  beams  A  and  B 
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into  four  beams,  which  are  nil  phase-shifted  interfere  grams  that  are  then  imaged  at  the 
sensing  camera.  The  opposite  handed  circularly  polarized  light  of  beam  A  develops  n 
shift  between  the  reflected  and  transmitted  beam  from  the  PBSl.  The  opposite  handed 
circularly  polarized  light  of  beam  B  develops  a  n  shift  between  the  reflected  and 
transmitted  beam  from  the  PBSl.  The  wave-front  phase  can  be  determined  from  the  four 
phase-shifted  images  using  a  four  bucket  method,  Equation  1  [22],  [23]. 

Figure  6  shows  an  example  of  the  four  beams  with  the  reference  only,  where  the 
signal  has  been  blocked,  and  the  signal  only,  where  the  reference  has  been  blocked. 
These  images  were  collected  with  a  16-bit  digital  imaging  system.  The  artificial  color 
mapping  from  blue  to  red  is  used  to  create  visual  distinction  between  low  and  high 
intensity  values  respectively. 


Reference  Only  Signal  Only 


Figure  6.  Image  of  four  phase  shifted  reference  beams,  and  four  phase  shifted  signal 
beams  before  interference.  Signal  beams  have  Rytov  number  of  0.044. 


This  SRI  utilizes  a  single  digital  imaging  system  to  collect  all  four  beams 
simultaneously  in  the  same  frame.  This  is  known  as  spatial  phase  shifting.  The  phase  of 
the  original  wave-front  can  then  be  extracted  from  the  image  of  these  four  beams  as 
previously  explained  [9]. 
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D. 


CAMERA 


The  camera  used  in  the  experimental  set-up  is  the  ALTA  U260  manufactured  by 
Apogee  Imaging  Systems.  This  camera  is  specifically  designed  for  low  noise  and  high 
dynamic  range  for  the  most  sensitive  of  scientific  imaging,  and  is  shown  in  Figure  7. 


Figure  7.  Apogee  Imaging  Systems  ALTA  U260  imaging  system  used  in  the  adaptive 

optics  test  bed  at  NFS. 

In  general,  the  signal  level  tends  to  drop  off  as  the  level  of  turbulence  increases. 
This  is  due  to  the  loss  of  coherence  in  the  beam,  and  the  dispersion  of  energy  as  the  beam 
propagates.  Integration  time  at  the  camera  is  selected  based  on  the  need  for  a  reasonable 
signal  level  at  the  highest  Rytov  numbers.  This  is  expected  to  yield  reasonable  SNR  for 
high  Rytov  numbers  and  more  than  sufficient  SNR  at  lower  Rytov  numbers.  In  a  closed- 
loop  adaptive  optics  setup  the  integration  time  would  be  a  small  fraction  of  this  value  due 
to  the  need  to  correct  the  wave-front  as  the  medium  itself  changes  at  a  high  frequency. 
The  details  of  the  imaging  system  are  shown  in  Table  2  [17]. 
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CCD 

Kodak  KAF- 
026  IE 

Digital  Resolution 

16  bits 

Array  Size  (pixels) 

512x512 

System 

Noise(typical) 

15  e'RMS 

Linear  Full  Well 
(typical) 

500,000  electrons 

Cooling(typical) 

50  C  below  ambient 

QE  at  632nm 

54% 

Dark 

Current(typical) 

1  eVpixel/see  (-25 

C) 

Table  2.  ALTA  U260  imaging  system  speeifieations. 
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IV.  RESULTS  AND  ANALYSIS 


The  experiment  was  eondueted  for  five  separate  Rytov  numbers  ranging  from 
0.044  to  0.4010.  This  range  was  neeessitated  by  the  experimental  set-up  eonstraints.  For 
Rytov  numbers  beyond  .4  the  signal  level  is  attenuated  too  mueh  to  get  useful  results. 
Each  experimental  result  will  be  shown  with  its  simulated  counterpart  for  the  purpose  of 
qualitative  analysis.  The  measure  of  effectiveness  of  the  experimental  results  will  be 
based  on  phase  reconstruction,  branch  point  identification,  and  overall  tilt.  This  data  is 
shown  in  Figures  8  to  22.  All  phase  data  is  reported  in  radians. 

A,  RYTOV  0.044 

1,  Phase  reconstruction 


Wrapped  Phase  Experimental  Wrapped  Phase  Simulated 

I 

■0 

I’ 

•D 

ma-2 


Figure  8.  Reconstructed  phase  data  for  experimental  and  simulated  cases  at  Rytov 

number  0.044. 


Mismatch  in  the  phase  reconstruction  can  be  explained  by  the  arbitrary  nature  of 
wrapping  cut  placement.  The  data  is  wrapped  between  -n  and  +n  due  to  the  domain  of 

the  tangent  function.  The  piston  term  in  the  wave-front  will  determine  how  far  forward 
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or  backward  the  wave-front  is  pushed,  and  determines  where  the  wrapping  cuts  lie.  In 
the  above  case  the  peak  of  the  simulated  data  has  been  wrapped,  but  in  the  experimental 
data  the  peak  remains  intact.  The  clear  consequence  of  this  is  that  the  experimental  data 
has  a  wrapping  cut  further  up  the  slope  toward  the  peak  of  the  data.  There  is  no  real 
impact  of  this  phenomenon  on  the  ability  of  the  adaptive  optics  system  to  compensate,  as 
piston  can  be  easily  removed  by  a  deformable. 

2,  Branch  Points 


Branch  Cuts  Experimental 


Branch  Cuts  Simulated 


Figure  9.  Branch  point  occurrence  for  experimental  and  simulated  cases  at  Rytov 

number  0.044. 


This  fairly  low  turbulence  level  did  not  generate  any  branch  points  in  the 
simulated  data.  The  occurrence  of  branch  points  in  the  experimental  data  around  the 
edges  can  be  explained  by  poor  quality  data  around  the  edges  of  the  beam  where  there 
may  be  incomplete  overlap  of  the  reference  and  signal  beams,  or  loss  of  signal  due  to 
optical  misalignment.  Branch  points  around  the  edges  of  the  data  are  a  common  theme 
throughout  the  experiment.  The  quality  of  this  reconstruction  is  quite  good.  There  are 
very  few  occurrences  of  branch  points  in  the  center  of  the  wave-front  where  the  data  is 
best. 
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3. 


Overall  Tilt 


Unwrapped  Phase  Experimental 


Unwrapped  Phase  Simulated 


Figure  10.  Unwrapped  phase  data  for  experimental  and  simulated  cases  at  Rytov 

number  0.044. 


After  the  phase  is  unwrapped  the  overall  tip/tilt  of  the  wave-front  can  be 
measured.  Piston  is  still  present,  as  in  the  wrapped  case,  but  the  important  factor  is  the 
total  range  of  the  data.  In  this  case  the  range  of  the  experimental  data  is  close  to  double 
that  of  the  simulated  data.  Once  again  the  edges  of  the  data  are  to  blame  because  they 
add  high  and  low  spikes  to  the  data  that  increase  the  overall  range.  If  these  outliers  are 
disregarded,  the  range  of  the  two  data  sets  is  quite  comparable. 
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B. 


RYTOV  0.109 


1,  Phase  Reconstruction 


Wrapped  Phase  Experimental 


Wrapped  Phase  Simulated 


Figure  1 1 .  Reconstructed  phase  data  for  experimental  and  simulated  cases  at  Rytov 

number  0.109. 

Increasing  the  turbulence  to  0.109  has  little  effect  qualitatively  on  wave-front 
reconstruction.  The  piston  term  still  shows  itself  in  the  variation  of  wrapping  cuts 
between  simulated  and  experimental  data.  The  overall  tilt  has  increased,  as  evidenced  by 
encroachment  of  the  wrapping  cut. 


24 


2. 


Branch  Points 


Branch  Cuts  Experimental  Branch  Cuts  Simulated 


Figure  12.  Branch  point  occurrence  for  experimental  and  simulated  cases  at  Rytov 

number  0.109. 

The  real  difference  that  can  be  seen  from  the  increased  level  of  turbulence  is  the 
addition  of  many  more  branch  points  in  the  experimental  data,  while  the  simulation  still 
shows  no  branch  cuts.  A  Rytov  number  of  0.109  is  not  yet  within  the  realm  of  deep 
turbulence.  As  turbulence  increases,  error  enters  in  the  form  of  misregistration,  signal 
drop-outs,  and  variation  in  alignment.  Misregistration  is  the  result  of  a  misalignment  of 
the  four  interfero grams,  which  must  be  compared  on  a  pixel  by  pixel  basis.  As 
turbulence  increases,  the  effects  of  misregistration  become  more  apparent  because  of  the 
increased  likelihood  that  adjacent  pixels  will  have  dramatically  different  phase  values, 
and  thus  will  skew  the  results  if  the  pixels  are  misaligned. 
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3. 


Overall  Tilt 


Unwrapped  Phase  Experimental 


--5 

-10 

1  I: 


Unwrapped  Phase  Simulated 


Figure  13.  Unwrapped  phase  data  for  experimental  and  simulated  eases  at  Rytov  number 

0.109. 


It  should  be  noted  that  the  overall  tilt  has  inereased  from  the  previous  Rytov 
number  from  approximately  8  to  around  12.  The  difference  in  overall  tilt  between  the 
experimental  and  simulated  data  is  still  explained  by  edge  effects,  though  with  greater 
difficulty.  The  peak  of  the  experimental  data  appears  to  be  quite  noisy,  and  the  spikes 
there  contribute  more  to  the  overall  tilt  than  edge  effects  can  explain. 
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c. 


RYTOV  0.205 


1,  Phase  Reconstruction 


Wrapped  Phase  Experimental  Wrapped  Phase  Simulated 

ft 


Figure  14.  Reconstructed  phase  data  for  experimental  and  simulated  cases  at  Rytov 

number  0.205. 

The  data  shows  an  additional  wrapping  cut  that  was  not  there  previously.  Also,  it 
is  becoming  apparent  that  the  experimental  data  has  a  decidedly  circular  appearance  as  if 
the  data  were  stretched  over  the  surface  of  a  sphere.  This  is  likely  due  to  the  camera 
being  out  of  the  imaging  plane  and  is  known  as  a  defocus  term.  This  is  the  second  order 
Zernike  polynomial  term.  Otherwise,  the  reconstruction  is  qualitatively  sound. 
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2. 


Branch  Points 


Branch  Cuts  Experimental 


Branch  Cuts  Simulated 


Figure  15.  Branch  point  occurrence  for  experimental  and  simulated  cases  at  Rytov 

number  0.205. 

This  Rytov  number  marks  the  appearance  of  a  branch  point  in  the  simulated  data. 
There  are  still  far  more  branch  points  in  the  experimental  data,  further  reinforcing  the 
suspicion  that  various  types  of  error  have  crept  in. 

3,  Overall  Tilt 


Unwrapped  Phase  Experimental 


Figure  16.  Unwrapped  phase  data  for  experimental  and  simulated  cases  at  Rytov  number 

0.205. 
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D. 


RYTOV  0.301 


1,  Phase  Reconstruction 


Wrapped  Phase  Experimental  Wrapped  Phase  Simulated 

I 

•0 

■ 

■0 

m^-2 


Figure  17.  Reconstructed  phase  data  for  experimental  and  simulated  cases  at  Rytov 

number  0.301. 


2,  Branch  Points 


Branch  Cuts  Experimental  Branch  Cuts  Simulated 


Figure  18.  Branch  point  occurrence  for  experimental  and  simulated  cases  at  Rytov 

number  0.301. 
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The  increased  level  of  turbulence  is  responsible  for  the  appearance  of  even  more 
branch  points  in  the  simulated  data.  The  turbulence  is  now  clearly  within  the  deep 
turbulence  regime.  The  number  of  branch  points  in  the  experimental  data  has  not 
significantly  increased,  which  may  indicate  the  experimental  data  is  showing  resilience  to 
increased  turbulence  levels.  There  are  for  the  first  significant  branch  cuts,  which  pose  a 
problem  for  traditional  reconstructors,  but  not  for  the  SRI. 

3,  Overall  Tilt 


Figure  19.  Unwrapped  phase  data  for  experimental  and  simulated  cases  at  Rytov  number 

0.301. 

The  experimental  and  simulated  data  are  remarkably  close  in  overall  tilt.  As  with 
the  branch  point  analysis,  this  may  indicate  the  experimental  set-up  is  showing  resilience 
to  increased  levels  of  turbulence. 


Unwrapped  Phase  Simulated 


Unwrapped  Phase  Experimental 
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E.  RYTOV  0.401 

1,  Phase  Reconstruction 

Wrapped  Phase  Experimental  Wrapped  Phase  Simulated 

I 

•0 

ig-2 

•0 

1^-2 

Figure  20.  Reconstructed  phase  data  for  experimental  and  simulated  cases  at  Rytov 

number  0.401. 

The  final  data  set  is  the  highest  level  of  turbulence  that  was  examined.  The 
number  of  wrapping  cuts  has  again  increased  indicating  an  increase  in  overall  tilt. 
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2. 


Branch  Points 


Branch  Cuts  Experimental 


Branch  Cuts  Simulated 


Figure  21 .  Branch  point  occurrence  for  experimental  and  simulated  cases  at  Rytov 

number  0.401. 


There  is  clearly  a  large  increase  in  branch  points  in  this  data  set.  The 
experimental  data  has  exploded  with  branch  points,  while  the  simulated  data  has  a  modest 
increase.  The  experimental  data  is  nearing  the  verge  of  its  capability  to  reproduce  the 
wave-front  accurately. 

3,  Overall  Tilt 


Unwrapped  Phase  Experimental 


Unwrapped  Phase  Simulated 


Figure  22.  Unwrapped  phase  data  for  experimental  and  simulated  cases  at  Rytov  number 

0.401. 
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V.  CONCLUSION 


The  SRI  was  successfully  developed  at  the  Naval  Postgraduate  School.  A 
qualitative  analysis  of  the  phase  reconstruction  in  the  experimental  results  shows  a 
remarkable  resemblance  to  the  phase  reconstruction  in  the  simulated  results.  Also,  upon 
comparison  of  the  unwrapped  phase  between  the  experimental  results  and  the  simulated 
results,  the  overall  tip  and  tilt  of  the  wave-front  is  found  to  be  comparable,  neglecting 
outliers.  Finally,  the  few  branch  points  that  do  appear  in  the  simulated  results  can  also  be 
identified  in  the  experimental  results.  The  similarity  of  these  features  between  the  ideal 
simulation  and  the  experimentation  demonstrates  the  success  of  the  experimental  set-up. 
While  there  is  still  significant  work  to  be  done  in  making  improvements  to  the  quality  of 
the  reconstruction  and  the  level  of  turbulence  the  system  can  handle,  the  instrument 
already  has  a  baseline  capability  to  reconstruct  wave-fronts  with  simulated  deep 
turbulence  effects.  Simulated  results  show  a  significant  drop-off  in  Strehl  ratio  below  a 
SNR  of  two.  The  system  is  limited  in  the  level  of  turbulence  that  can  be  accurately 
processed  due  to  loss  of  signal  strength  at  the  sensor.  Beyond  a  Rytov  number  of  0.4,  the 
system  creeps  out  of  alignment  and  the  reconstruction  can  no  longer  be  completed 
because  of  diminishing  SNR. 

Further  improvements  need  to  be  made  to  the  instrument  before  it  can  be  used  in 
closed-loop  adaptive  optics  and  beam  correction.  Further  characterization  should  be  done 
on  the  spatial  light  modulator,  and  an  effort  made  to  identify  error  sources  in  the 
simulation  of  deep  turbulence.  The  system  as  a  whole  should  have  its  alignment  verified 
and  brought  into  focus  to  eliminate  Fresnel  ringing  in  the  reference  path,  and  to  bring  the 
camera  into  the  focal  plane  of  the  interferogram.  Consideration  should  be  given  to 
replacing  the  spatial  filter  with  a  single-mode  optical  amplifier  for  amplification  of  the 
reference  path.  This  should  also  reduce  the  systems  sensitivity  to  overall  tip  and  tilt  in 
the  signal,  which  alters  the  alignment  of  the  beam  and  causes  it  to  “miss”  the  spatial  filter 
on  the  reference  path  of  the  SRI. 

Future  work  includes  using  closing  the  adaptive  optics  loop  with  a  deformable 

mirror  to  apply  beam  correction,  incorporating  a  Shack-Hartmann  WFS  to  compliment 
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the  SRI,  testing  the  SRI  with  higher  levels  of  turbulence,  and  developing  control 
algorithms  for  the  deformable  mirror  that  are  specifically  designed  for  SRI  phase 
reconstruction  data. 
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APPENDIX 


A,  THESIS  MAIN  SCRIPT 

close  all; 
clear  all; 
clc; 


%DATA  FORMATTING 

g, 

0 

a 

0 

%Experimental  Data 

%Pass  along  the  image  you  wish  to  parse  from  the  workspace: 
load ( ' S ; \MATLAB\Thesis\Experimental_Data_Final .mat ' ) 

Rytov  =  0.401; 

[Fl_Exp, F2_Exp, F3_Exp, F4_Exp, dim]  =  Manual_Acquire (Rytov_0401_Both) ; 

g, 

0 

%Simulated  Data 

%Pass  along  Rytov  number  of  requested  data  from  following  list: 

%  0.044;  0.109;  0.205;  0.301;  0.401;  0.610;  1.090;  1.308;  1.613;  2.092; 
%  2.310;  2.615;  3.051;  3.313;  3.618;  4.010 

[Fl_Sim, F2_Sim, F3_Sim, F4_Sim, AMP_Sim, PHA_Sim]  =  SRI_Sim (Rytov, dim) ; 

q, 

0 

%Camera  Noise 
SNR  =  4; 

[Fl_Sim,Nl]  =  Noise (le-3, Fl_Sim, SNR) ; 

[F2_Sim,N2]  =  Noise ( le-3 , F2_Sim, SNR) ; 

[F3_Sim,N3]  =  Noise ( le-3 , F3_Sim, SNR) ; 

[F4_Sim,N4]  =  Noise ( le-3 , F4_Sim, SNR) ; 

%Plot_Frames (N1 , N2 , N3 ,  N4 ) 

q, 

0 

%Generate  Frame  Plots 

Plot_Frames (Fl_Exp, F2  Exp,F3  Exp,F4  Exp) 

Plot_Frames (Fl_Sim,  F2_Sim, F3_Sim, F4_Sim) 

g, 

0 

%Four  Bucket  Method 

Phase  Wrapped_Exp  =  atan2 ( (F4  Exp  -  F2  Exp)  ,  (Fl_Exp  -  F3  Exp) ) ; 

Phase  Wrapped_Sim  =  atan2 ( (F4  Sim  -  F2_Sim) , (Fl_Sim  -  F3_Sim) ) ; 

Complex  Wave_Exp  =  AMP_Sim. *exp ( j . *Phase_Wrapped_Exp) ; 

Complex  Wave_Sim  =  AMP  Sim. *exp ( j . *Phase  Wrapped_Sim) ; 

q, 

0 

q, 

o 


%STREHL  RATIO  COMPUTATION 

q, 

0 

q, 

0 

Uorig  =  AMP  Sim. *exp ( j . *PHA  Sim) ; 

Urec_Sim  =  exp (- j . *Phase_Wrapped_Sim) ; 

%how  well  I  reconstructed  the  phase,  not  as  meaningful  for  real  data 
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%because  I  don't  have  the  original  phase/amp  for  this  data  set.  Always 
one 

%if  no  noise  added. 

Strehl_Reconstructor_Sim  = 

(abs  ( sum  ( sum  (Uorig  .  *Urec_Sim)  ))  ''2)  /  { sum  ( sum  (abs  (Uorig)  )  )  )  ''2; 

q, 

0 

o 

0 


%UNWRAP 

q, 

0 

g, 

o 

%Unwrap  Mask 

%Select  circular  mask  or  square  mask  for  unwrapping 
%For  square  mask  set  mask_flag  =  0; 

%For  circular  mask  set  mask_flag  =  1; 
mask_flag  =  1; 

if  (mask_flag  ==  0) 

Mask  =  ones (dim) ; 

end 

if  (mask_flag  ==  1) 

Mask  =  zeros (dim) ; 
for  i=  1 : dim 

for  j  =  1 : dim 

r  =  norm ( [round (dim/2 )  round (dim/2 )]-[ i  j]); 
if  (r  <  dim/2) 

Ma  s  k ( i , j )  =  1  ; 

end 

end 

end 

end 

g, 

0 

%Unwrap  Algorithm 

%Unwrap  the  Phase  using  one  of  the  below  algorithms 
%For  Quality  Guided  2D  Unwrap  method  =  1; 

%For  Goldstein  2D  Unwrap  method  =2; 

%For  Constantin!  Network  Programming  2D  Unwrap  method  =  3; 
method  =2; 

Phase  Unwrapped_Sim  =  Phase_Unwrap (  Complex  Wave_Sim,  Mask,  method); 
Phase_Unwrapped_Exp  =  Phase_Unwrap (Complex_Wave_Exp,  Mask,  method); 

q, 

0 

q, 

o 


B.  FRAME  ACQUISITION  FUNCTION 

function  [FI , F2 , F3 , F4 , dim]  =  Manual_Acquire ( Im_in) 

%MANUAL_ACQUIRE  Takes  in  the  full  camera  image  and  return  4  sub  frames. 
%The  frames  have  been  adjusted  for  beam  size,  beam  shape, 

%and  average  intensity. 
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%specify  the  image  you  wish  to  process 
Image  =  Im_in; 


%Begin  Frame  Grabbing 

%Optimized  for  New  SLM  image,  must  re-optimize  for  new  alignment 
F3  =  Image ( 35 : 242 , 270 : 4 67 )  ; 

F4  =  Image ( 50 : 2 60 , 56 : 254 )  ; 

F2  =  Image (2 82 : 4 95 , 73 : 2 68 ) ; 

FI  =  Image (262 : 475, 295 : 500)  ; 
dim  =  200; 


FI 

=  imresize (FI , 

[dim 

dim] ) ; 

F2 

=  imresize (F2 , 

[dim 

dim] ) ; 

F3 

=  imresize (F3, 

[dim 

dim] ) ; 

F4 

=  imresize (F4 , 

[dim 

dim] ) ; 

IM 

mask  =  zeros (dim) ; 

for 

i=  1 : dim 

for  j  =  1 : dim 

r  =  norm ( [round (dim/2 )  round (dim/2 )]- [i  j ] )  ; 
if  (r  <  dim/2) 

IM_ma  s  k ( i , j )  =  1 ; 

end 

end 

end 

%normalize  the  average  intensity  of  the  four  frames. 
Fl_temp  =  F1.*IM  mask; 

F2  temp  =  F2.*IM  mask; 

F3_temp  =  F3.*IM_mask; 

F4_temp  =  F4.*IM  mask; 


Fl_sum 
F2  sum 
F3_sum 
F4  sum 


sum (sum (Fl_temp) ) ; 
sum (sum (F2_temp) ) ; 
sum (sum (F3_temp) ) ; 
sum (sum (F4_temp) ) ; 


Fl_avg 

F2_avg 

F3_avg 

F4_avg 


Fl_sum  / (pi* (dim/2 ) ^2 ) ; 
F2_sum  /  (pi*  (dim/2  )  ^^2  )  ; 
F3_sum  / (pi* (dim/2 ) ^2 ) ; 
F4  sum  / (pi* (dim/2) ^2) ; 


Avg_avg  =  (Fl_avg  +  F2_avg  +  F3_avg  +  F4_avg)/4; 


FI  =  FI . *Avg_avg . /Fl_avg; 
F2  =  F2 . *Avg_avg . /F2  avg; 
F3  =  F3 . *Avg_avg . /F3  avg; 
F4  =  F4 . *Avg_avg . /F4_avg; 

end 
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c 


SRI  SIMULATOR  FUNCTION 


function  [  FI , F2 , F3 , F4 , AMP, PHA  ]  =  SRI_Sim(  Rytov,dim  ) 
%SRI_SIM  Summary  of  this  function  goes  here 
%  Detailed  explanation  goes  here 

load ( ' S ; \MATLAB\Thesis\Simulation  Data  Final. mat') 


%The  following  data  is  used  in  the  simulation  and  experimentation 
portion 

%of  the  Experiment 

Q-Q-9-9-9-Q-9-Q-9-9-9-9-9-9-9-Q-9-9-9-9-9-9-9-Q-9-9-9-9-9-9-9-Q-9-9-9-9-9-9-9-Q-9-9-9-Q-9-9-9-Q-9-9-9-9-9-9-9-Q-9-9-9-9-9-9-9-9-9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

9-  9-  9-  9- 
0  0  0  0 

%Rytov  0.044 
if  (Rytov  ==  0.044) 

FI  =  AMP_0044  +  AMP_0044 . *cos (PHA_0044) ; 

F2  =  AMP_0044  -  AMP_0044 . *sin (PHA_0044) ; 

F3  =  AMP_0044  -  AMP_0044 . *cos (PHA_0044) ; 

F4  =  AMP_0044  +  AMP_0044 . *sin (PHA_0044) ; 

AMP  =  AMP_0044; 

PHA  =  PHA_0044; 

end 


%Rytov  0.109 
if  (Rytov  ==  0.109) 

FI  =  AMP_0109  +  AMP_0109 . *cos (PHA_0109)  ; 
F2  =  AMP_0109  -  AMP_0109 . *sin (PHA_0109) ; 
F3  =  AMP_0109  -  AMP_0109 . *cos (PHA_0109) ; 
F4  =  AMP_0109  +  AMP_0109 . *sin (PHA_0109) ; 
AMP  =  AMP_0109; 

PHA  =  PHA_0109; 

end 


%Rytov  0.205 
if  (Rytov  ==  0.205) 

FI  =  AMP_0205  +  AMP_0205 . *cos (PHA_0205)  ; 
F2  =  AMP_0205  -  AMP_0205 . *sin (PHA_0205) ; 
F3  =  AMP_0205  -  AMP_0205 . *cos (PHA_0205) ; 
F4  =  AMP_0205  +  AMP_0205 . *sin (PHA_0205) ; 
AMP  =  AMP_0205; 

PHA  =  PHA_0205; 

end 


%Rytov  0.301 
if  (Rytov  ==  0.301) 

FI  =  AMP_0301  +  AMP_0301 . *cos (PHA_0301)  ; 
F2  =  AMP_0301  -  AMP_0301 . *sin (PHA_0301) ; 
F3  =  AMP_0301  -  AMP_0301 . *cos (PHA_0301) ; 
F4  =  AMP_0301  +  AMP_0301 . *sin (PHA_0301) ; 
AMP  =  AMP_0301; 

PHA  =  PHA_0301; 

end 
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o\o  o\o 


%Rytov  0.401 
if  (Rytov  ==  0.401) 


FI  = 

AMP 

0401  + 

AMP 

0401 . *cos (PHA 

0401) ; 

F2  = 

AMP 

0401  - 

AMP 

’0401 .  *sin  (PHA 

'0401)  ; 

F3  = 

AMP 

0401  - 

AMP 

0401 . *cos (PHA 

^0401)  ; 

F4  = 

AMP 

0401  + 

AMP 

’0401 .  *sin  (PHA 

^0401)  ; 

AMP 

=  AMP 

0401; 

PHA 

=  PHA 

0401; 

end 

9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9- 

00000000000000000000000000000000000000000000000000000000000000000000000 

g,  g,  g,  g, 
o  o  0  o 


The  following  data  is  used  only  in  the  Simulation  portion  of  the 
Experiment 

%Rytov  0.610 
if  (Rytov  ==  0.610) 


FI  = 

AMP 

0610  + 

AMP 

0610 . *cos (PHA 

0610 

F2  = 

AMP 

0610  - 

AMP 

0610  .  *sin (PHA 

^0610 

F3  = 

AMP 

0610  - 

AMP 

0610 . *cos (PHA 

'0610 

F4  = 

AMP 

0610  + 

AMP 

0610 . *sin (PHA 

’0610 

AMP 

=  AMP 

0610; 

PHA 

=  PHA 

0610; 

end 


%Rytov  1.090 
if  (Rytov  ==  1.090) 


FI  = 

AMP 

1090  + 

AMP 

1090 . *cos 

(PHA 

1090 

F2  = 

AMP 

1090  - 

AMP 

1090 . *sin 

(PHA 

'l090 

F3  = 

AMP 

1090  - 

AMP 

1090 . *cos 

(PHA 

1090 

F4  = 

AMP 

1090  + 

AMP 

1090  .  *sin 

(PHA 

^1090 

AMP 

=  AMP 

1090; 

PHA 

=  PHA 

1090; 

end 


%Rytov  1.308 
if  (Rytov  ==  1.308) 
FI  =  AMP_1308  + 
F2  =  AMP_1308  - 
F3  =  AMP_1308  - 
F4  =  AMP_1308  + 
AMP  =  AMP_1308; 
PHA  =  PHA_1308; 

end 


AMP_1308 . *cos (PHA_1308) ; 
AMP_1308 . *sin (PHA_1308) ; 
AMP_1308 . *cos (PHA_1308) ; 
AMP  1308 . *sin (PHA  1308); 


%Rytov  1.613 
if  (Rytov  ==  1.613) 


FI  = 

AMP 

1613  + 

AMP 

1613 . *cos 

(PHA 

1613 

F2  = 

AMP 

1613  - 

AMP 

’l613.*sin 

(PHA 

'l613 

F3  = 

AMP 

1613  - 

AMP 

1613 . *cos 

(PHA 

^1613 

F4  = 

AMP 

1613  + 

AMP 

1613 . *sin 

(PHA 

^1613 

AMP 

=  AMP 

1613; 

PHA 

=  PHA 

1613; 

end 
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%Rytov  2.092 


(Rytov  == 

2.092) 

FI  =  AMP 

2092  + 

AMP 

2092 . *cos 

(PHA 

F2  =  AMP 

'2092  - 

AMP 

2092 . *sin 

(PHA 

F3  =  AMP 

'2092  - 

AMP 

2092 . *cos 

(PHA 

F4  =  AMP 

'2092  + 

AMP 

2092  .  *sin 

(PHA 

AMP  =  AMP_2092; 
PHA  =  PHA_2092; 

end 

%Rytov  2.310 
if  (Rytov  ==  2.310) 


FI 

=  AMP 

2310 

+ 

AMP 

2 

3 

10 . *cos 

(PHA 

F2 

=  AMP 

'2310 

- 

AMP 

2 

3 

10  .  *sin 

(PHA 

F3 

=  AMP 

'2310 

- 

AMP 

'2 

3 

10 . *cos 

(PHA 

F4 

=  AMP 

'2310 

+ 

AMP 

'2 

3 

10  .  *sin 

(PHA 

AMP  =  AMP_2310; 
PHA  =  PHA_2310; 

end 

%Rytov  2.615 
if  (Rytov  ==  2.615) 


FI 

=  AMP 

2615 

+ 

AMP 

2615 . *cos 

(PHA 

F2 

=  AMP 

'2615 

- 

AMP 

2615 . *sin 

(PHA 

F3 

=  AMP 

'2615 

- 

AMP 

2615 . *cos 

(PHA 

F4 

=  AMP 

'2615 

+ 

AMP 

2615 . *sin 

(PHA 

AMP  =  AMP_2615; 
PHA  =  PHA_2615; 

end 

%Rytov  3.051 
if  (Rytov  ==  3.051) 


FI 

=  AMP 

3051 

+ 

AMP 

3051 . *cos 

(PHA 

F2 

=  AMP 

'3051 

- 

AMP 

3051 . *sin 

(PHA 

F3 

=  AMP 

'3051 

- 

AMP 

3051 . *cos 

(PHA 

F4 

=  AMP 

'3051 

+ 

AMP 

3051 . *sin 

(PHA 

AMP  =  AMP_3051; 
PHA  =  PHA  3051; 


end 


%Rytov  3.313 

if  (Rytov  == 

3.313) 

FI  =  AMP 

3313  + 

AMP 

3313 . *cos 

(PHA 

F2  =  AMP 

3313  - 

AMP 

3313 . *sin 

(PHA 

F3  =  AMP 

3313  - 

AMP 

3313 . *cos 

(PHA 

F4  =  AMP 

3313  + 

AMP 

'3313. *sin 

(PHA 

AMP  =  AMP 

‘  3313; 

PHA  =  PHA 

._3313; 

end 

%Rytov  3.618 

if  (Rytov  == 

3.618) 

FI  =  AMP 

3618  + 

AMP 

3618 . *cos 

(PHA 

F2  =  AMP 

3618  - 

AMP 

3618 . *sin 

(PHA 

2092) ; 
2092) ; 
2092) ; 
2092) ; 


2310) ; 
2310) ; 
2310) ; 
2310) ; 


2615) ; 
2615) ; 
2615) ; 
2615) ; 


3051) ; 
3051)  ; 
3051)  ; 
3051) ; 


3313) ; 
3313) ; 
3313)  ; 
3313) ; 


3618)  ; 
3618)  ; 
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F3  =  AMP_3618  -  AMP_3618 . *cos (PHA_3618) ; 
F4  =  AMP_3618  +  AMP_3618 . *sin (PHA_3618) ; 
AMP  =  AMP_3618; 

PHA  =  PHA_3618; 

end 


%Rytov  4.010 
if  (Rytov  ==  4.010) 


FI  = 

AMP 

4010  + 

AMP 

4010 . *cos (PHA 

4010)  ; 

F2  = 

AMP 

4010  - 

AMP 

4010  .  *sin (PHA 

*4010)  ; 

F3  = 

AMP 

4010  - 

AMP 

4010  .  *cos (PHA 

’4010)  ; 

F4  = 

AMP 

4010  + 

AMP 

4010 . *sin (PHA 

’4010)  ; 

AMP 

=  AMP 

4010; 

PHA 

=  PHA 

4010; 

end 


%Resize  and  mirror  data  to  more  closely  match  experimental 
FI  =  imresize(Fl,  [dim  dim]); 

F2  =  imresize(F2,  [dim  dim]); 

F3  =  imresize(F3,  [dim  dim]); 

F4  =  imresize(F4,  [dim  dim]); 

AMP  =  imresize (AMP,  [dim  dim]); 

[dim  dim] ) ; 


data 


PHA  = 
Tempi 
Temp2 
Temp  3 
Temp  4 
Temp  5 
Temp  6 
for  i 


imresize (PHA, 
=  zeros (dim) ; 
=  zeros (dim) ; 
=  zeros (dim) ; 
=  zeros (dim) ; 
=  zeros (dim) ; 
=  zeros (dim) ; 
=  1 : dim 

,  i) 

,  i) 

,  i) 

,  i) 

,  i) 

,  i) 


Tempi ( 
Temp2 ( 
Temp3 ( 
Temp4 ( 
TempS ( 
Temp6 ( 


FI  ( 
F2  ( 
F3( 
F4  ( 


, abs ( i- 
, abs ( i- 
, abs ( i- 
, abs ( i- 


(dim+1 ) 
(dim+1 ) 
(dim+1 ) 
(dim+1 ) 


AMP ( : , abs ( i- (dim+1 ) ) ) ; 
PHA ( : , abs ( i- (dim+1 ) ) ) ; 


end 

FI  =  Tempi; 
F2  =  Temp2; 
F3  =  Temp3; 
F4  =  Temp4; 
AMP  =  TempS; 
PHA  =  Temp 6; 


end 


D.  NOISE  FUNCTION 

function  [  Noisy_Image,  Noise_Frame]  =  Noise (  Int_time, Image, SNR) 

%NOISE  :  Noise  modeled  after  the  ALTA  U260  camera  from  Apogee  Imaging 
%  Systems,  which  contains  the  KAF-0261  CCD  from  TrueSense  imaging,  inc 
%  This  imaging  system  is  a  cooled  system,  with  low  noise,  and  high 
dynamic 
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%  range . 

%Definition  of  noise  in  this  case  is  SNR  =  mu/sigma,  which  is  the 
%mean  of  the  signal  divided  by  the  std  dev  of  the  noise. 

%Define  Camera  Parameters 
[rows,-]  =  size ( Image) ; 

signal_average  =  sum ( sum ( Image) ) /rows ^2 ; 

sigma_total_noise  =  signal_average/SNR;  %this  includes  all  system  noise 
dark_current_rate  =  10;  %this  is  for  -25  C 

non_unif ormity  =  .05;  %standard  deviation  is  5%  of  pixel  value 
[rows,-]  =  size (Image); 

dark_current_f rame  =  zeros (rows) +  dark  current_rate*Int_time; 

bias_frame  =  sigma_total_noise . *randn (rows) ; 

non_uniform  =  1  +  non  uniformity . *randn (rows) ; 

Noise_Frame  =  dark_current_f rame  +  bias_frame; 

Noisy_Image  =  Image . *non_uniform  +  Noise_Frame; 

E.  2D  PHASE  UNWRAP  FUNCTION 


9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9- 

00000000000000000000000000000000000000000000000000000000000000000000000 

Q-Q-9-Q-9-Q-9-9-9-9-9-9-9-Q-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-Q-9- 

ooooooooooooooooooooooooooooooo 

%  GoldsteinUnwrap2D  is  a  script  to  demonstrate  the  2D  Goldstein  branch 
cut  phase  unwrapping  algorithm. 

g, 

o 

%  Calls:  PhaseResidues 
%  BranchCuts 

%  FloodFill 

g, 

0 

%  References : : 

%  1.  R.  M.  Goldstein,  H.  A.  Zebken,  and  C.  L.  Werner,  "Satellite  radar 
interferometry : 

%  Two-dimensional  phase  unwrapping,"  Radio  Sci.,  vol.  23,  no.  4,  pp . 
713-720,  1988. 

%  2.  D.  C.  Ghiglia  and  M.  D.  Pritt,  Two-Dimensional  Phase  Unwrapping: 

%  Theory,  Algorithms  and  Software.  New  York:  Wiley-Interscience, 

1998  . 

g, 

0 

%  Inputs:  1.  IM  =  Complex  image 
%  Optional  inputs: 

%  2.  im_mask  =  Binary  mask 

%  3.  max_box_radius=  Maximum  search  box  radius  (pixels) 

%  4.  threshold  std  =  Number  of  noise  standard  deviations  used 


for 


%  thresholding  the  magnitude  image 

%  Outputs:  1.  Unwrapped  phase  image 
%  2.  Phase  quality  map 

g, 

0 

%  This  code  can  easily  be  extended  for  3D  phase  unwrapping. 

g, 

o 

%  Posted  by  Bruce  Spottiswoode  on  22  December  2008 
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%  2010/07/23  Modified  by  Carey  Smith 

%  1.  Changed  IM_mask,  IM_mag,  IM  phase  to  lowercase,  in 

order  to  be 

%  similar  to  the  Quality  Guided  routines. 

%  2 .  Moved  the  logic  to  chose  the  reference  point  from 

FloodFill .m 

%  to  here,  in  order  to  be  similar  to  the  Quality  Guided 

routines . 

%  The  coordinates  of  the  reference  point,  colref, 

rowref,  are  now 

%  passed  to  FloodFill2.m 

%  3.  User  can  optionally  have  the  code  automatically  chose 

the 


%  largest  magnatude  point  as  the  reference  point. 

%  4.  Allowed  the  user  to  specify  IM,  im  mask, 

max_box_radius  before 

%  calling  this  routine.  (Note  that  threshold_std  is  not 


used . ) 


%  5.  Modified  the  plots  to  suit  my  preferences. 

g, 

o 
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0000000000000000000000000000000 


IM  =  Complex_Wave; 

im  mag  =  abs(IM);  %Magnitude  image 

im_phase  =  angle (IM);  %Phase  image 

%%  Replace  with  your  mask  (if  desired) 

im  mask  =  Mask; 

mag_max  =  max ( im  mag ( : ) ) ; 

%indxl  =  find(im_mag  <  0.1*mag  max);  %Intensity  =  mag^2,  so  this  =  .04 
threshold  on  the  intensity 

%im  mask(indxl)  =  0;  %  Don't  mask  at  this  point;  wait  until  residues 

are  computed 

if ( ~exist ( ' im  mask ' , ' var ' ) ) 

im_mask  =  ones (size (IM) ) ;  %Mask  (if  applicable) 

end 

figure;  imagesc (im_mag. *im  mask),  colormap (gray) ,  axis  square,  axis 
off,  title ('GS  Initial  masked  magnitude');  colorbar; 

figure;  imagesc (im_phase . *im_mask) ,  colormap (gray) ,  axis  square,  axis 
off,  title ('GS  Initial  masked  phase');  colorbar; 

%%  Compute  the  residues 

residue_charge  =  PhaseResidues  rl (im_phase,  im_mask) ;  %  Calculate  phase 
residues  (Does  not  use  mask) 

figure;  imagesc (residue_charge) ,  colormap (gray) ,  axis  square,  axis  off, 
title ('GS  Phase  residues  (charged) ');  colorbar; 

%%  Compute  the  branch  cuts 

max_box_radius=floor (length (residue_charge) /2) ;  %  Maximum  search  box 

radius  (pixels) 

%max_box  radius=4  %  Maximum  search  box  radius  (pixels) 
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if ( ~exist ( ' max_box_radius ' , ' var ' ) ) 

max_box_radius=4 ;  %  Maximum  search  box  radius  (pixels) 

end 

%  BranchCutsO  ignores  residues  with  mask  ==  0,  so  keep  the  entire  mask 

branch_cuts  =  BranchCuts_rl (residue_charge,  max_box_radius ,  im_mask) ;  % 
Place  branch  cuts 

figure;  imagesc (branch_cuts ) ,  colormap (gray) ,  axis  square,  axis  off, 
title ('GS  Branch  cuts');  colorbar; 

im  mask (branch_cuts )  =  0;  %  Now  need  to  mask  off  branch  cut  points,  in 

order  to  avoid  an  error  in  FloodFill 

im_magl  =  im_mag . *im_mask;  %  Mask  off  magnitude  ==  0  points,  so  that 
they  are  not  chosen  for  the  starting  point 

%%  Manually  (default)  or  automatically  identify  starting  seed  point 
if(l)  %  Chose  starting  point  interactively 
im_phase_quality  =  im_magl; 

minp  =  im_phase_quality (2 : end-1 ,  2:end-l);  minp  =  min (minp ( : ) ) ; 
maxp  =  im  phase_quality (2 : end-1 ,  2:end-l);  maxp  =  max (maxp ( : ) ) ; 
figure;  imagesc ( im_phase_quality,  [minp  maxp]),  colormap (gray)  , 
colorbar,  axis  square,  axis  off;  title (' Phase  quality  map'); 

%uiwait (msgbox (' Select  known  true  phase  reference  phase  point.  Black 
=  high  quality  phase;  white  =  low  quality  phase.', 'Phase  reference 
point ' , ' modal ' ) ) ; 

uiwait (msgbox (' Select  known  true  phase  reference  phase  point.  White  = 
high  magnitude;  Black  =  low  magnitude.', 'Phase  reference 
point ' , 'modal ' ) ) ; 

[xpoint, ypoint]  =  ginput(l);  %Select  starting  point  for  the 

guided  floodfill  algorithm 
colref  =  round (xpoint) ; 
rowref  =  round (ypoint) ; 

close;  %Close  the  figure; 

else  %  Chose  starting  point  =  max.  intensity 
[r_dim,  c_dim] =size (im  phase); 

im  magl(l, :)  =0;  %Set  magnitude  of  border  pixels 

to  0,  so  that  they  are  not  used  for  the  reference 
im_magl (r_dim, : )  =0; 
im_magl ( : , 1 )  =  0 ; 
im  magi ( : , c_dim)  =  0; 

[rowrefn, colrefn]  =  find(im  magi  >=  0.99*mag  max); 

rowref  =  rowrefn (1);  %Choose  the  1st  point  for  a 

reference  (known  good  value) 

colref  =  colrefn(l);  %Choose  the  1st  point  for  a 

reference  (known  good  value) 
end 


%%  Unwrap 

if (exist ( ' rowref '  ,  ' var ' ) ) 

im_unwrapped  =  FloodFill_rl (im_phase,  im_mag,  branch_cuts,  im_mask, 
colref,  rowref) ;  %  Flood  fill  phase  unwrapping 
else 

im_unwrapped  =  FloodFill_rl (im_phase,  im_mag,  branch_cuts,  im_mask) ; 
%  Flood  fill  phase  unwrapping 
end 
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axis  square. 


%  Display  results 

figure;  imagesc ( im_unwrapped) ,  colormap (gray) , 
axis  off,  title ('GS  Unwrapped  phase'); 


colorbar , 
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